{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Microsoft Differential Privacy Whitepaper Collateral Notebooks Part 3\n",
    "# Privacy-Preserving Statistical Analysis with Differential Privacy "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/histogram.jpg\" width=800 />\n",
    "\n",
    "The goal of this notebook is to demonstrate how to generate and release basic statistical outcomes in a differntially private fashion. We will assess the impact on various privacy levels and dataset sizes on accuracy.\n",
    "\n",
    "We are using a dataset conataining demographic data about Californian residents from the Public Use Microdata Sample (PUMS) statistics.\n",
    "The dataset includes more than 1.2m records and is therefore well suited to experiment with different sample sizes.\n",
    "\n",
    "The data includes personal information like gender, age, ethnical background, income and marital status. For our analysis, we will focus on creating histograms for the yearly income.\n",
    "\n",
    "The default method for generating a histogram in SmartNoise is by releasing counts of each bin or category using the geometric mechanism. The geometric mechanism only returns integer values for any query, so resists some vulnerabilities of DP releases from floating point approximations. It is also possible to generate histograms from the more typical Laplace mechanism. We show both approaches below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installs and imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %pip install opendp seaborn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "import seaborn as sns\n",
    "\n",
    "datafile = os.path.join('.','data','PUMS-california-demographics.csv')\n",
    "samplefile = os.path.join('.','data','sample.csv') # export generated sample and import for DP analysis\n",
    "\n",
    "df = pd.read_csv(datafile, usecols = range(1,11))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>state</th>\n",
       "      <th>puma</th>\n",
       "      <th>sex</th>\n",
       "      <th>age</th>\n",
       "      <th>educ</th>\n",
       "      <th>income</th>\n",
       "      <th>latino</th>\n",
       "      <th>black</th>\n",
       "      <th>asian</th>\n",
       "      <th>married</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>6</td>\n",
       "      <td>60100</td>\n",
       "      <td>0</td>\n",
       "      <td>83</td>\n",
       "      <td>9</td>\n",
       "      <td>20500.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>6</td>\n",
       "      <td>60100</td>\n",
       "      <td>1</td>\n",
       "      <td>81</td>\n",
       "      <td>9</td>\n",
       "      <td>4800.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>6</td>\n",
       "      <td>60100</td>\n",
       "      <td>0</td>\n",
       "      <td>45</td>\n",
       "      <td>9</td>\n",
       "      <td>12000.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6</td>\n",
       "      <td>60100</td>\n",
       "      <td>1</td>\n",
       "      <td>42</td>\n",
       "      <td>12</td>\n",
       "      <td>7200.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>6</td>\n",
       "      <td>60100</td>\n",
       "      <td>0</td>\n",
       "      <td>35</td>\n",
       "      <td>11</td>\n",
       "      <td>55600.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>6</td>\n",
       "      <td>60100</td>\n",
       "      <td>1</td>\n",
       "      <td>34</td>\n",
       "      <td>11</td>\n",
       "      <td>4600.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>6</td>\n",
       "      <td>60100</td>\n",
       "      <td>0</td>\n",
       "      <td>20</td>\n",
       "      <td>8</td>\n",
       "      <td>8000.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>6</td>\n",
       "      <td>60100</td>\n",
       "      <td>0</td>\n",
       "      <td>35</td>\n",
       "      <td>12</td>\n",
       "      <td>43000.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   state   puma  sex  age  educ   income  latino  black  asian  married\n",
       "0      6  60100    0   83     9  20500.0       0      0      0        1\n",
       "1      6  60100    1   81     9   4800.0       0      0      0        1\n",
       "2      6  60100    0   45     9  12000.0       0      0      0        1\n",
       "3      6  60100    1   42    12   7200.0       0      0      0        1\n",
       "4      6  60100    0   35    11  55600.0       0      0      0        1\n",
       "5      6  60100    1   34    11   4600.0       0      0      0        1\n",
       "6      6  60100    0   20     8   8000.0       1      0      0        0\n",
       "7      6  60100    0   35    12  43000.0       0      0      0        0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocess data\n",
    "Converting the continuous variable income to catagories as foundation for the further analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove entries with negative income values\n",
    "df.drop (df[df['income'] < 0].index, inplace=True)\n",
    "\n",
    "income_categories = ['less 10k', '10k-20k', '20k-30k', '30k-40k', '40k-50k', '50k-60k', '60k-80k', '80k-100k', '100k-150k', 'above 150k']\n",
    "cut_bins = [-1, 10000, 20000, 30000, 40000, 50000, 60000, 80000, 100000, 150000, 1000000] # 10 categories\n",
    "df['inc_cat'] = pd.cut(df['income'], bins=cut_bins, labels=income_categories)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Review histogram of income from original dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABdoAAAPeCAYAAAAbM6gSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAABYlAAAWJQFJUiTwAAByqklEQVR4nOz9fZRtdX3n+36+QguKCYrtSXwIByGiJtE2G7kqRvDhajQhIVE84DkdiURt7YAJoOkcRMU0enICwajYci9pJWrfAV68wYPtw01EBd0GZGtC0qRBga2tMUFFQR4T5Hv+WLMyVpa1N7X3r4qi4PUaY42511zzO+eviv/eNZmrujsAAAAAAMDOud96LwAAAAAAADYyoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCO/9CVf2Xqvov670OAAAAAICNYtf1XgD3OI/btGnTpiT/63ovBAAAAADgblY7M+SOdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMCAXdd7Adw7HPC69633EtgBW0596XovAQAAAADuNdzRDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBgw4f2qvq3VdXT6+XbOObQqvp0Vd1QVTdV1SVVddRdnPeoqrp0Ov6Gaf7Q7Ry/S1UdV1WXV9WtVXV9VX20qg7azswDqurNVXVlVd1WVddV1Qer6vHbmdmrqv6oqrZW1e1V9XdV9Z6qetT2fh4AAAAAANbGhg7tVfUTSc5IctN2jjkmyQVJfibJB5KcleQRSc6uqtO2MXNakrOTPHw6/gNJnpDkgul8i8dXknOSnJ7k/tOa/jTJwUkuqqrDlpnZLcmfJXljkhuTvD3Jnyf51SSXVdVTlpl5aJLPJ/mtJFcneVuSS5O8LMmWqtp3W78HAAAAAADWxq7rvYCdNcXt9yb5TpL/X5LXLnPMPklOS3J9kid399Zp/+8l+UKSE6rqQ939+bmZg5KckFnIPrC7vzvtPzXJliSnVdVHls41OTLJ4Uk2J3lOd982zZyZ5LNJzqqqC7v7+3Mzxyd5epLzkhzR3XdOM+cmOT/Je6rqCUv7J29Nsn+S07v7hLk1vyazUP+fkjx/Jb8/AAAAAABWx0a+o/01SZ6d2d3cN2/jmKOT7JbkjPkwPsXzt05vX7Uws/T+LUuRfZrZmuRd0/letjDz6ml70lJkn2a+kOTcJA/LLMQn+ec/Eixd53fmY3p3fzjJxUl+KskhczMPSvJr08968sL1z0jy1SQ/7652AAAAAIC714YM7dMzzH8/ydu7+6LtHPrsafvxZT772MIxOzVTVbsnOSjJLZkF8pVcZ78keye5qruvXeHMU5M8IMnnFu6MzxTqPzG9fdYy5wMAAAAAYI1suEfHVNWuSd6f5GtJTryLwx87ba9a/KC7v1lVNyd5VFU9sLtvqao9kjwyyU3d/c1lzvflabv/3L79kuyS5JruvmOFM9tc1yrPbFNVbdnGR49byTwAAAAAADMbLrRn9uWhP5vk57r71rs4ds9pe8M2Pr8hyR7Tcbes8PgkefAOXmO9ZgAAAAAAWGMbKrRX1VMyu4v9D+e/wJQd190HLLd/utN90928HAAAAACADWvDPKN9emTM+zJ7dMobVji2dJf3ntv4fPEu8ZUe/72duMZ6zAAAAAAAsMY2TGhP8qDMnj/++CS3VVUvvZK8aTrmrGnfH03vr5y2P/Tc8qp6eGaPjfl6d9+SJN19c5JvJHnQ9Pmix0zb+eekX53kB0n2nf4YsJKZba5rlWcAAAAAAFhjG+nRMbcn+c/b+GxTZs9t/2xmQXrpsTIXJnl6kufP7Vvygrlj5l2Y5Nemmffe1Ux331ZVm5M8Y3p9agXXuTqzL3Pdv6oe3d3XrmDmL5LcmuTpVfUj3f39pQ+q6n5Jnje9Xbw+AAAAAABraMPc0d7dt3b3y5d7Jfm/psP+ZNp37vT+vZkF+mOqap+lc1XVQzJ71nuSnLlwqaX3r5+OW5rZJ8lvTudbDPDvnranVNXuczMHJjkiybeSfGjuZ+m56/zBFMqXZg7LLNhfkeQzczM3JXl/Znfhn7xw/WOS7JPkE919TQAAAAAAuNtspDvad1h3X1tVr0vyjiSXVdW5Sf4xyeFJHpVlvlS1uzdX1elJjk9yeVWdl+T+mQXzvZIc291bFy51TpIXTuf9UlVdkOSh08wuSV7R3TcuzJye5NBp5pKq+mSSvZO8OMktSY7u7jsXZk5M8swkx1fVk5JcmtmjdA5Lcl1mfwgAAAAAAOButGHuaN9Z3f3OJL+c5L8leWmSVyb5+yS/3t2v3cbMCUleNh33ymnuvyX5pe4+Y5njO8lLMovzdyQ5NrPwflGSg7v7w8vM3J7kuUn+Y5IHJzluen9+kgO7+5JlZr6T5GmZ/eHgJ5OckOQpmd1hf0B3X72CXwkAAAAAAKuoZo0YZqpqy6ZNmzZt2bJlh+YOeN371mhFrIUtp750vZcAAAAAAPdEtTND9/o72gEAAAAAYC0J7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABiw4UJ7Vf2fVfXJqvofVXVrVV1fVV+qqjdV1UMXjt2nqno7r3O2c52jqurSqrqpqm6oqk9X1aHbOX6Xqjquqi6fW9dHq+qg7cw8oKreXFVXVtVtVXVdVX2wqh6/nZm9quqPqmprVd1eVX9XVe+pqkfd1e8OAAAAAIDVt+t6L2AnHJfki0n+LMl1SfZI8tQkJyd5ZVU9tbv/x8LMXyU5f5lz/c1yF6iq05KckOTrSc5Kcv8kRya5oKqO7e4zFo6vJOckOTzJlUnOSLJXkiOSXFRVL+ruDy/M7Db9DE9PclmStyf5iSQvTvKLVfXs7r5kYeahSTYn2T/JhdM1H5fkZdPM07r7muV+JgAAAAAA1sZGDO0/2t23Le6sqrckOTHJ/57k3y98/JfdffJKTj7dgX5CkquTHNjd3532n5pkS5LTquoj3b11buzIzCL75iTPWVpfVZ2Z5LNJzqqqC7v7+3Mzx2cW2c9LckR33znNnJvZHwXeU1VPWNo/eWtmkf307j5hbs2vySzU/6ckz1/JzwkAAAAAwOrYcI+OWS6yTz44bR8zeIlXTdu3LEX26bpbk7wryW6Z3UE+79XT9qT59XX3F5Kcm+RhmYX4JP98B/zSdX5nPqZPd75fnOSnkhwyN/OgJL+W5ObM7t6fd0aSryb5+arad+U/KgAAAAAAozZcaN+OX5q2ly/z2SOq6t9V1YnT9onbOc+zp+3Hl/nsYwvHpKp2T3JQklsyC+R3OZNkvyR7J7mqu69d4cxTkzwgyecW7ozPFOo/Mb191jLnAwAAAABgjWzER8ckSarqtUkelGTPJE9O8nOZRfbfX+bw506v+flPJzmqu782t2+PJI9MclN3f3OZ83x52u4/t2+/JLskuaa771jhzGOn7VXLHL+aM9tUVVu28dHjVjIPAAAAAMDMhg3tSV6b5Mfm3n88ya9397fm9t2S5D9m9szzpS8JfWJmj155VpJPVtWTuvvm6bM9p+0N27jm0v4Hz+27J88AAAAAALDGNmxo7+4fT5Kq+rHMHt3y+0m+VFWHdvcXp2OuS/LGhdGLqup5mX1J6VOSvDyzLxK9T+nuA5bbP93pvuluXg4AAAAAwIa14Z/R3t3/0N1/muR5SR6a5H0rmLkjyR9Pbw+e+2jprvA9s7yl/d/bIDMAAAAAAKyxDR/al3T3V5NckeSnq+pfr2Bk6REze8yd4+Yk30jyoKp6+DIzj5m2889JvzrJD5LsW1XL/R8Cy81cOW239Tz11ZoBAAAAAGCN3WtC++QR0/YHKzj2qdP2moX9F07b5y8z84KFY9LdtyXZnOSBSZ6xkpnM4vzXkuxfVY9e4cxfJLk1ydOr6kfmD66q+2V2R3+SfGqZ8wEAAAAAsEY2VGivqv2r6ocenVJV96uqtyT5n5Js7u7vTvs3TRF68fjnJDluevuBhY/PnLavr6qHzM3sk+Q3k9ye5L0LM++etqdU1e5zMwcmOSKzu+c/tLS/u3vuOn8wv8aqOiyzYH9Fks/MzdyU5P2Z3YF/8sL1j0myT5JPdPfiHw4AAAAAAFhDG+3LUH8hyf9RVZ9Ncm2S7yT5sSSHJNk3yd8necXc8acneUxVbU7y9WnfE5M8e/r3G7p78/wFuntzVZ2e5Pgkl1fVeUnun1kw3yvJsd29dWFd5yR5YZLDM/tC1gsye178EUl2SfKK7r5xYeb0JIdOM5dU1SeT7J3kxUluSXJ0d9+5MHNikmcmOb6qnpTk0iSPT3JYkusy+0MAAAAAAAB3o40W2v88yU8m+bkkP5vkwUluzuy55O9P8o7uvn7u+Pcn+dUkB2b2OJZ/leQfknwwyRndffFyF+nuE6rqrzML169McmeSLyY5tbs/sszxXVUvyewRMkcnOTbJbUkuSnLKYsyfZm6vqucm+d0kL8nsDvsbk5yf5E3dfcUyM9+pqqcleVOSX8nszvfvZHaH/Ru7++uLMwAAAAAArK2aPcUEZqpqy6ZNmzZt2bJlh+YOeN371mhFrIUtp750vZcAAAAAAPdEtTNDG+oZ7QAAAAAAcE8jtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGDAhgvtVfV/VtUnq+p/VNWtVXV9VX2pqt5UVQ/dxsxBVfXR6dhbq+ryqvrtqtplO9c5tKo+XVU3VNVNVXVJVR11F2s7qqounY6/YZo/dDvH71JVx03rWfpZPlpVB21n5gFV9eaqurKqbquq66rqg1X1+O2tDQAAAACAtbHhQnuS45LskeTPkrw9yX9JckeSk5NcXlU/MX9wVR2W5KIkByf50yRnJLl/krclOWe5C1TVMUkuSPIzST6Q5Kwkj0hydlWdto2Z05KcneTh0/EfSPKEJBdM51s8vqbrnz6t54xpfQcnuWha9+LMbtPP/cYkN04//58n+dUkl1XVU5ZbGwAAAAAAa6e6e73XsEOqavfuvm2Z/W9JcmKSd3f3v5/2/WiSryTZM8nTu/uypXMkuTDJ05K8pLvPmTvPPkn+e5KbkxzQ3Vun/Q9J8oUk+yU5qLs/PzdzUJLPJbk6yYHd/d25c23J7A8Dj1s61/TZS5L8f5JsTvKcpZ+pqg5M8tkkNyTZr7u/Pzfzvyd5a5LzkhzR3XdO+w9Lcn6SK5I8YWn/zqiqLZs2bdq0ZcuWHZo74HXv29lLsg62nPrS9V4CAAAAANwT1c4Mbbg72peL7JMPTtvHzO07PMnDkpyzFNnnznHS9PbVC+c5OsluSc6YD+NTPH/r9PZVCzNL79+yFNmnma1J3jWd72ULM0vXPWn+Z+ruLyQ5d1r34Uv7pzvgl67zO/Mxvbs/nOTiJD+V5JAAAAAAAHC32XChfTt+adpePrfv2dP248scf1GSW5IcND2SZSUzH1s4ZqdmpjvqD5quf/EKr7Nfkr2TXNXd1+7A2pZVVVuWeyV53ErmAQAAAACY2XW9F7Czquq1SR6U2WNhnpzk5zKL7L8/d9hjp+1Vi/PdfUdVXZvkp5Psm+RvVzDzzaq6OcmjquqB3X1LVe2R5JFJburuby6z1C9P2/3n9u2XZJck13T3HSuc2ea6tjMDAAAAAMAa27ChPclrk/zY3PuPJ/n17v7W3L49p+0N2zjH0v4H7+DMHtNxt6zhNVZjZpu6+4Dl9k93tW9ayTkAAAAAANjAj47p7h/v7kry40lemNld6V+qKpEYAAAAAIC7zYYN7Uu6+x+6+0+TPC/JQ5O8b+7jpbu89/yhwX+5/3s7MXPDwnYtrjE6AwAAAADAGtvwoX1Jd381yRVJfrqq/vW0+8pp+0PPLa+qXZM8OskdSa6Z+2h7Mw/P7LExX+/uW6br3pzkG0keNH2+6DHTdv7Z6lcn+UGSfad1rGRmm+vazgwAAAAAAGvsXhPaJ4+Ytj+YthdO2+cvc+zBSR6YZHN33z63f3szL1g4Zqdmuvu2JJun6z9jhde5OsnXkuxfVY/egbUBAAAAALCGNlRor6r9q+qHHp1SVferqrck+Z8yC+ffnT46L8m3kxxZVU+eO373JKdMb9+9cLr3Jrk9yTFVtc/czEOSnDi9PXNhZun966fjlmb2SfKb0/neuzCzdN1TpvUszRyY5Igk30ryoaX93d1z1/mDqrrf3MxhmQX7K5J8JgAAAAAA3G2We2zJPdkvJPk/quqzSa5N8p0kP5bkkMy+DPXvk7xi6eDuvrGqXpFZcP90VZ2T5Pokv5zksdP+c+cv0N3XVtXrkrwjyWVVdW6Sf0xyeJJHJfnD7v78wszmqjo9yfFJLq+q85LcP7NgvleSY7t768LPck5mX+J6eGZf4npBZs+YPyLJLkle0d03LsycnuTQaeaSqvpkkr2TvDjJLUmO7u47V/KLBAAAAABgdWy00P7nSX4yyc8l+dkkD05yc2bPJX9/knd09/XzA919flUdkuT1SV6UZPckX8ksir9julM8CzPvrKqtSV6b5KWZ3fl/RZKTuvtPlltYd59QVX+d2R3sr0xyZ5IvJjm1uz+yzPFdVS/J7BEyRyc5NsltSS5Kckp3b15m5vaqem6S303ykiTHJbkxyflJ3tTdVyz7WwMAAAAAYM3UMp2Z+7Cq2rJp06ZNW7Zs2aG5A173vjVaEWthy6kvXe8lAAAAAMA9Ue3M0IZ6RjsAAAAAANzTCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAAzZUaK+qh1bVy6vqT6vqK1V1a1XdUFWfrarfqKr7LRy/T1X1dl7nbOdaR1XVpVV103SNT1fVods5fpeqOq6qLp/WdX1VfbSqDtrOzAOq6s1VdWVV3VZV11XVB6vq8duZ2auq/qiqtlbV7VX1d1X1nqp61F39/gAAAAAAWH27rvcCdtCLk7w7yTeTfCrJ15L8WJIXJvnjJC+oqhd3dy/M/VWS85c5398sd5GqOi3JCUm+nuSsJPdPcmSSC6rq2O4+Y+H4SnJOksOTXJnkjCR7JTkiyUVV9aLu/vDCzG5J/izJ05NcluTtSX5i+hl/saqe3d2XLMw8NMnmJPsnuXC65uOSvGyaeVp3X7PczwQAAAAAwNrYaKH9qiS/nOS/dvedSzur6sQklyZ5UWbR/UMLc3/Z3Sev5ALTHegnJLk6yYHd/d1p/6lJtiQ5rao+0t1b58aOzCyyb07ynO6+bZo5M8lnk5xVVRd29/fnZo7PLLKfl+SIpZ+nqs7N7I8C76mqJ8z/nEnemllkP727T5hb82syC/X/KcnzV/JzAgAAAACwOjbUo2O6+8LuvmAhPqe7/z7JmdPbZw5e5lXT9i1LkX26xtYk70qyW2Z3kM979bQ9aSmyTzNfSHJukodlFuKT/PMd8EvX+Z35n2e68/3iJD+V5JC5mQcl+bUkNyc5eeH6ZyT5apKfr6p9V/6jAgAAAAAwakOF9rvwT9P2jmU+e0RV/buqOnHaPnE753n2tP34Mp99bOGYVNXuSQ5KcktmgfwuZ5Lsl2TvJFd197UrnHlqkgck+dzCnfGZQv0nprfPWuZ8AAAAAACskY326JhlVdWuSV46vV0ukD93es3PfDrJUd39tbl9eyR5ZJKbuvuby5zny9N2/7l9+yXZJck13b1c5F9u5rHT9qpljl/NmW2qqi3b+OhxK5kHAAAAAGDm3nJH++8n+ZkkH+3uT8ztvyXJf0xyQJKHTK9DMvsi1Wcm+eQU15fsOW1v2MZ1lvY/eIPMAAAAAACwxjb8He3TF4GekOS/Z/YM83/W3dcleePCyEVV9bzMvqT0KUlentkXid6ndPcBy+2f7nTfdDcvBwAAAABgw9rQd7RX1TGZRfIrkjyru69fydz0iJc/nt4ePPfR0l3he2Z5S/u/t0FmAAAAAABYYxs2tFfVbyd5Z5K/ySyy//0OnuJb0/afHx3T3Tcn+UaSB1XVw5eZecy0nX9O+tVJfpBk3+lZ8SuZuXLabut56qs1AwAAAADAGtuQob2q/kOStyX5y8wi+3U7cZqnTttrFvZfOG2fv8zMCxaOSXfflmRzkgcmecZKZjKL819Lsn9VPXqFM3+R5NYkT6+qH5k/uKrul+R509tPLXM+AAAAAADWyIYL7VX1hsy+/HRLkud097e3c+ymKUIv7n9OkuOmtx9Y+PjMafv6qnrI3Mw+SX4zye1J3rsw8+5pe0pV7T43c2CSIzK7e/5DS/u7u+eu8wfza6yqwzIL9lck+czczE1J3p/ZHfgnL1z/mCT7JPlEdy/+4QAAAAAAgDW0ob4MtaqOSvJ7mT2q5eIkr6mqxcO2dvfZ079PT/KYqtqc5OvTvicmefb07zd09+b54e7eXFWnJzk+yeVVdV6S+2cWzPdKcmx3b1245jlJXpjk8CRfqqoLkjx0mtklySu6+8aFmdOTHDrNXFJVn0yyd5IXJ7klydHdfefCzIlJnpnk+Kp6UpJLkzw+yWFJrsvsDwEAAAAAANyNNlRoT7L0mJVdkvz2No75TJKzp3+/P8mvJjkws8ex/Ksk/5Dkg0nO6O6LlztBd59QVX+dWbh+ZZI7k3wxyand/ZFlju+qeklmj5A5OsmxSW5LclGSUxZj/jRze1U9N8nvJnlJZnfY35jk/CRv6u4rlpn5TlU9LcmbkvxKZne+fyezO+zf2N1fX5wBAAAAAGBt1ewpJjBTVVs2bdq0acuWLTs0d8Dr3rdGK2ItbDn1peu9BAAAAAC4J/qhR6isxIZ7RjsAAAAAANyTCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAgFUN7VW1d1X96F0c8yNVtfdqXhcAAAAAANbLat/Rfm2S37qLY14zHQcAAAAAABveaof2ml4AAAAAAHCfsB7PaP/xJDevw3UBAAAAAGDV7Tp6gqp66cKuJy2zL0l2SbJ3kn+b5K9HrwsAAAAAAPcEw6E9ydlJevp3Jzlsei1aeqTMLUnevArXBQAAAACAdbcaof1l07aSvCfJ+Uk+vMxxP0jynSSf7+7vrcJ1AQAAAABg3Q2H9u7+k6V/V9VRSc7v7veNnhcAAAAAADaC1bij/Z9197NW83zAxve133vCei+BHbD3G32FBgAAAMCOut96LwAAAAAAADayVQ/tVXVIVX2kqq6rqn+qqh8s87pjta8LAAAAAADrYVUfHVNVv5jZl6HukuRrSa5MIqoDAAAAAHCvtaqhPcnJSf4pyS929/9/lc8NAAAAAAD3OKv96JifSXKuyA4AAAAAwH3Faof2m5Jcv8rnBAAAAACAe6zVDu2fTPK0VT4nAAAAAADcY612aP8PSfarqpOqqlb53AAAAAAAcI+z2l+G+qYk/y3Jm5McXVV/meR7yxzX3f0bq3xtAAAAAAC42612aP/1uX/vM72W00mEdgAAAAAANrzVDu2PXuXzAQAAAADAPdqqhvbu/upqng8AAAAAAO7pVvvLUAEAAAAA4D5lVe9or6q9V3psd39tNa8NAAAAAADrYbWf0b41sy86vSu9BtcGAAAAAIC73WrH7vdl+dD+4CRPSvI/J/l0Es9yBwAAAADgXmG1vwz117f1WVXdL8kbkrwqyVGreV0AAAAAAFgvd9uXoXb3nd395sweL/P7d9d1AQAAAABgLd1toX3O5iTPW4frAgAAAADAqluP0L5Xkj3W4boAAAAAALDq7tbQXlX/zyRHJPmbu/O6AAAAAACwVlb1y1Cr6sLtXOcnkuw9vf+91bwuAAAAAACsl1UN7UmeuY39neS7ST6R5LTu3laQBwAAAACADWVVQ3t3r8cz3wEAAAAAYN0I4wAAAAAAMGC1Hx3zL1TVjyR5cJIbuvvGtbwWAAAAAACsh1W/o72qdq2q362qryT5XpKtSb5bVV+Z9q9p3AcAAAAAgLvTqkbvqrp/ko8nOSSzL0D9H0m+meThSfZJ8pYkz6+q53X3P67mtQEAAAAAYD2s9h3txyd5ZpL/muTx3b1Pdz+tu/dJ8tgkFyR5xnQcAAAAAABseKsd2v/XJH+T5Fe6+8vzH3T31UlemOS/JfnfVvm6AAAAAACwLlY7tP9kko91953LfTjt/1iS/Vb5ugAAAAAAsC5WO7T/Y5IH3cUxeyT5p1W+LgAAAAAArIvVDu2XJzm8qh623IdV9a+THJ7kr1b5ugAAAAAAsC5WO7SfkeRhSS6tqt+oqn2r6gFV9eiqelmSS6bPz1jl6wIAAAAAwLrYdTVP1t0frKonJfndJP/vZQ6pJH/Q3R9czesCAAAAAMB6WdXQniTdfWJV/V9JfiPJzybZM8kNSb6U5D3d/fnVviYAAAAAAKyXVQ/tSdLdf5HkL9bi3AAAAAAAcE+yqs9or6oXV9WFVfWIbXz+yKr6ZFW9cDWvCwAAAAAA62W1vwz15Uke3N1/t9yH3f2NzB4l8/JVvi4AAAAAAKyL1Q7tT0hy2V0c84UkT1zl6wIAAAAAwLpY7dC+V5Lr7uKY7yT516t8XQAAAAAAWBerHdq/neQxd3HMY5J8b5WvCwAAAAAA62K1Q/vnkvxyVT1uuQ+r6vFJDkty8SpfFwAAAAAA1sVqh/bTkuya5LNV9Zqq2r+q9pi2v5VZYN9lOg4AAAAAADa8XVfzZN39har690neleRt02veD5K8ursvWc3rAgAAAADAelnV0J4k3X1WVX02yb9P8pQkD87smex/keTd3f23q31NAAAAAABYL6se2pNkiunHrsW5AQAAAADgnmS1n9EOAAAAAAD3KUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGbKjQXlUPraqXV9WfVtVXqurWqrqhqj5bVb9RVcv+PFV1UFV9tKqun2Yur6rfrqpdtnOtQ6vq09P5b6qqS6rqqLtY31FVdel0/A3T/KHbOX6XqjpuWs+t0/o+WlUHbWfmAVX15qq6sqpuq6rrquqDVfX47a0NAAAAAIC1saFCe5IXJzkryVOSXJLkj5J8KMnPJPnjJB+sqpofqKrDklyU5OAkf5rkjCT3T/K2JOcsd5GqOibJBdN5PzBd8xFJzq6q07Yxc1qSs5M8fDr+A0mekOSC6XyLx9d0/dOn9Zwxre/gJBdN616c2S3JnyV5Y5Ibk7w9yZ8n+dUkl1XVU5ZbGwAAAAAAa2fX9V7ADroqyS8n+a/dfefSzqo6McmlSV6U5IWZxfdU1Y9mFr1/kOSZ3X3ZtP8NSS5McnhVHdnd58yda58kpyW5PsmTu3vrtP/3knwhyQlV9aHu/vzczEFJTkhydZIDu/u70/5Tk2xJclpVfWTpXJMjkxyeZHOS53T3bdPMmUk+m+Ssqrqwu78/N3N8kqcnOS/JEUu/g6o6N8n5Sd5TVU+Y/90AAAAAALC2NtQd7d19YXdfsBiSu/vvk5w5vX3m3EeHJ3lYknOWIvt0/G1JTprevnrhMkcn2S3JGfNhfIrnb53evmphZun9W5Yi+zSzNcm7pvO9bGFm6bonLUX2aeYLSc6d1n340v7pDvil6/zO/O+guz+c5OIkP5XkkAAAAAAAcLfZUKH9LvzTtL1jbt+zp+3Hlzn+oiS3JDloeiTLSmY+tnDMTs1U1e5JDpquf/EKr7Nfkr2TXNXd1+7A2gAAAAAAWEMb7dExy6qqXZO8dHo7H7sfO22vWpzp7juq6tokP51k3yR/u4KZb1bVzUkeVVUP7O5bqmqPJI9MclN3f3OZ5X152u4/t2+/JLskuaa77/jhkWVntrmu7cxsU1Vt2cZHj1vJPAAAAAAAM/eWO9p/P7MvLv1od39ibv+e0/aGbcwt7X/wTszsubBdi2uMzgAAAAAAsMY2/B3tVfWazL6I9L8n+bV1Xs6G0d0HLLd/utN90928HAAAAACADWtD39FeVcckeXuSK5I8q7uvXzhk8e7zRUv7v7cTMzcsbNfiGqMzAAAAAACssQ0b2qvqt5O8M8nfZBbZ/36Zw66ctj/03PLpue6PzuzLU69Z4czDk+yR5OvdfUuSdPfNSb6R5EHT54seM23nn61+dZIfJNl3WsdKZra5ru3MAAAAAACwxjZkaK+q/5DkbUn+MrPIft02Dr1w2j5/mc8OTvLAJJu7+/YVzrxg4Zidmunu25Jsnq7/jBVe5+okX0uyf1U9egfWBgAAAADAGtpwob2q3pDZl59uSfKc7v72dg4/L8m3kxxZVU+eO8fuSU6Z3r57Yea9SW5PckxV7TM385AkJ05vz1yYWXr/+um4pZl9kvzmdL73LswsXfeUaT1LMwcmOSLJt5J8aGl/d/fcdf6gqu43N3NYZsH+iiSfCQAAAAAAd5sN9WWoVXVUkt/L7LErFyd5TVUtHra1u89Oku6+sapekVlw/3RVnZPk+iS/nOSx0/5z54e7+9qqel2SdyS5rKrOTfKPSQ5P8qgkf9jdn1+Y2VxVpyc5PsnlVXVekvtnFsz3SnJsd29dWOc5SV44nfdLVXVBkodOM7skeUV337gwc3qSQ6eZS6rqk0n2TvLiJLckObq779zuLxEAAAAAgFW1oUJ7Zs9UT2Yh+re3ccxnkpy99Ka7z6+qQ5K8PsmLkuye5CuZRfF3THeK/wvd/c6q2prktUlemtmd/1ckOam7/2S5i3b3CVX115ndwf7KJHcm+WKSU7v7I8sc31X1ksweIXN0kmOT3JbkoiSndPfmZWZur6rnJvndJC9JclySG5Ocn+RN3X3FNn4nAAAAAACskQ0V2rv75CQn78Tc55L8wg7OXJDkgh2cOTtzkX8Fx9+R2bPm37YDM7ckeeP0AgAAAABgnW24Z7QDAAAAAMA9idAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAgF3XewEA3Dc9/Z1PX+8lsAM+d+zn1nsJAAAAcI/ljnYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAgA0X2qvq8Kp6Z1VdXFU3VlVX1Qe2cew+0+fbep2znescVVWXVtVNVXVDVX26qg7dzvG7VNVxVXV5Vd1aVddX1Uer6qDtzDygqt5cVVdW1W1VdV1VfbCqHr+dmb2q6o+qamtV3V5Vf1dV76mqR21rBgAAAACAtbPrei9gJ5yU5N8kuSnJ15M8bgUzf5Xk/GX2/81yB1fVaUlOmM5/VpL7JzkyyQVVdWx3n7FwfCU5J8nhSa5MckaSvZIckeSiqnpRd394YWa3JH+W5OlJLkvy9iQ/keTFSX6xqp7d3ZcszDw0yeYk+ye5cLrm45K8bJp5Wndfs4LfBwAAAAAAq2QjhvbjMgvgX0lySJJPrWDmL7v75JWcfLoD/YQkVyc5sLu/O+0/NcmWJKdV1Ue6e+vc2JGZRfbNSZ7T3bdNM2cm+WySs6rqwu7+/tzM8ZlF9vOSHNHdd04z52b2R4H3VNUTlvZP3ppZZD+9u0+YW/NrMgv1/ynJ81fycwIAAAAAsDo23KNjuvtT3f3l7u41usSrpu1bliL7dN2tSd6VZLfM7iCf9+ppe9JSZJ9mvpDk3CQPyyzEJ/nnO+CXrvM78zF9uvP94iQ/ldkfEpZmHpTk15LcnOTkheufkeSrSX6+qvZd+Y8KAAAAAMCoDRfad9IjqurfVdWJ0/aJ2zn22dP248t89rGFY1JVuyc5KMktmQXyu5xJsl+SvZNc1d3XrnDmqUkekORzC3fGZwr1n5jePmuZ8/2Qqtqy3CsrexQPAAAAAACTjfjomJ3x3On1z6rq00mO6u6vze3bI8kjk9zU3d9c5jxfnrb7z+3bL8kuSa7p7jtWOPPYaXvVNta7WjMAAAAAAKyxe3tovyXJf8zsmedLXxL6xMwevfKsJJ+sqid1983TZ3tO2xu2cb6l/Q+e23dPntmm7j5guf3TXe2bVnIOAAAAAADu5Y+O6e7ruvuN3f3F7v7e9LooyfOSXJLkJ5O8fH1XCQAAAADARnavDu3bMj3i5Y+ntwfPfbR0V/ieWd7S/u9tkBkAAAAAANbYfTK0T741bfdY2jE9QuYbSR5UVQ9fZuYx03b+OelXJ/lBkn2rarlH8Sw3c+W03dbz1FdrBgAAAACANXZfDu1PnbbXLOy/cNo+f5mZFywck+6+LcnmJA9M8oyVzGQW57+WZP+qevQKZ/4iya1Jnl5VPzJ/cFXdL7PH4STJp5Y5HwAAAAAAa+ReHdqratMUoRf3PyfJcdPbDyx8fOa0fX1VPWRuZp8kv5nk9iTvXZh597Q9pap2n5s5MMkRmd09/6Gl/d3dc9f5g/k1VtVhmQX7K5J8Zm7mpiTvz+wO/JMXrn9Mkn2SfKK7F/9wAAAAAADAGlruUSf3aFX1K0l+ZXr749P2aVV19vTvb3f3a6d/n57kMVW1OcnXp31PTPLs6d9v6O7N8+fv7s1VdXqS45NcXlXnJbl/ZsF8ryTHdvfWhWWdk+SFSQ5P8qWquiDJQ6eZXZK8ortvXJg5Pcmh08wlVfXJJHsneXGSW5Ic3d13LsycmOSZSY6vqicluTTJ45McluS6zP4QAAAAAADA3WjDhfYkT0py1MK+fadXknw1yVJof3+SX01yYGaPY/lXSf4hyQeTnNHdFy93ge4+oar+OrNw/cokdyb5YpJTu/sjyxzfVfWSzB4hc3SSY5PcluSiJKcsxvxp5vaqem6S303ykszusL8xyflJ3tTdVywz852qelqSN2X2x4ZnJPlOZnfYv7G7v744AwAAAADA2tpwob27T84PPzplW8f+5yT/eSevc3aSs3fg+DuSvG16rXTmliRvnF4rnbk+yW9NLwAAAAAA1tm9+hntAAAAAACw1oR2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAbsut4LAACY95mDD1nvJbCDDrnoM+u9BAAAgHXljnYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYsOFCe1UdXlXvrKqLq+rGquqq+sBdzBxUVR+tquur6taquryqfruqdtnOzKFV9emquqGqbqqqS6rqqLu4zlFVdel0/A3T/KHbOX6XqjpuWs+t0/o+WlUHbWfmAVX15qq6sqpuq6rrquqDVfX47a0NAAAAAIC1seFCe5KTkhyT5ElJvnFXB1fVYUkuSnJwkj9NckaS+yd5W5JztjFzTJILkvxMkg8kOSvJI5KcXVWnbWPmtCRnJ3n4dPwHkjwhyQXT+RaPr+n6p0/rOWNa38FJLprWvTizW5I/S/LGJDcmeXuSP0/yq0kuq6qn3NXvAwAAAACA1bXrei9gJxyX5OtJvpLkkCSf2taBVfWjmUXvHyR5ZndfNu1/Q5ILkxxeVUd29zlzM/skOS3J9Ume3N1bp/2/l+QLSU6oqg919+fnZg5KckKSq5Mc2N3fnfafmmRLktOq6iNL55ocmeTwJJuTPKe7b5tmzkzy2SRnVdWF3f39uZnjkzw9yXlJjujuO6eZc5Ocn+Q9VfWEpf0AAAAAAKy9DXdHe3d/qru/3N29gsMPT/KwJOcsRfbpHLdldmd8krx6YeboJLslOWM+jE/x/K3T21ctzCy9f8tSZJ9mtiZ513S+ly3MLF33pKXIPs18Icm507oPX9o/3QG/dJ3fmY/p3f3hJBcn+anM/vgAAAAAAMDdZMOF9h307Gn78WU+uyjJLUkOmh7JspKZjy0cs1MzVbV7koOm61+8wuvsl2TvJFd197U7sDYAAAAAANbQRnx0zI547LS9avGD7r6jqq5N8tNJ9k3ytyuY+WZV3ZzkUVX1wO6+par2SPLIJDd19zeXWcOXp+3+c/v2S7JLkmu6+44VzmxzXduZ2aaq2rKNjx63knkAAAAAAGbu7Xe07zltb9jG50v7H7wTM3subNfiGqMzAAAAAACssXv7He1sQ3cfsNz+6U73TXfzcgAAAAAANqx7+x3ti3efL1ra/72dmLlhYbsW1xidAQAAAABgjd3bQ/uV0/aHnlteVbsmeXSSO5Jcs8KZhyfZI8nXu/uWJOnum5N8I8mDps8XPWbazj9b/eokP0iy77SOlcxsc13bmQEAAAAAYI3d20P7hdP2+ct8dnCSBybZ3N23r3DmBQvH7NRMd9+WZPN0/Wes8DpXJ/lakv2r6tE7sDYAAAAAANbQvT20n5fk20mOrKonL+2sqt2TnDK9fffCzHuT3J7kmKraZ27mIUlOnN6euTCz9P7103FLM/sk+c3pfO9dmFm67inTepZmDkxyRJJvJfnQ0v7u7rnr/EFV3W9u5rDMgv0VST4TAAAAAADuNhvuy1Cr6leS/Mr09sen7dOq6uzp39/u7tcmSXffWFWvyCy4f7qqzklyfZJfTvLYaf+58+fv7mur6nVJ3pHksqo6N8k/Jjk8yaOS/GF3f35hZnNVnZ7k+CSXV9V5Se6fWTDfK8mx3b114Uc5J8kLp/N+qaouSPLQaWaXJK/o7hsXZk5Pcug0c0lVfTLJ3klenOSWJEd3953b/QUCAAAAALCqNlxoT/KkJEct7Nt3eiXJV5O8dumD7j6/qg5J8vokL0qye5KvZBbF3zHdKf4vdPc7q2rrdJ6XZnbn/xVJTuruP1luUd19QlX9dWZ3sL8yyZ1Jvpjk1O7+yDLHd1W9JLNHyByd5NgktyW5KMkp3b15mZnbq+q5SX43yUuSHJfkxiTnJ3lTd1+x3NoAAAAAAFg7Gy60d/fJSU7ewZnPJfmFHZy5IMkFOzhzdpKzd+D4O5K8bXqtdOaWJG+cXgAAAAAArLN7+zPaAQAAAABgTQntAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwYNf1XgAAAKzUGSdcsN5LYAcc84e/tN5LAACAu4U72gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADNh1vRcAAAAAa+Vv33Lhei+BHfD41z97vZcAADvFHe0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYsOt6LwAAAGDUW/7t4eu9BHbA6z9w3novAQBgVbmjHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAy4T4T2qtpaVb2N199vY+agqvpoVV1fVbdW1eVV9dtVtct2rnNoVX26qm6oqpuq6pKqOuou1nZUVV06HX/DNH/odo7fpaqOm9Zz67S+j1bVQSv/jQAAAAAAsFp2Xe8F3I1uSPJHy+y/aXFHVR2W5ENJbktybpLrk/xSkrcleXqSFy8zc0ySdyb5TpIPJPnHJIcnObuqntDdr11m5rQkJyT5epKzktw/yZFJLqiqY7v7jIXjK8k503mvTHJGkr2SHJHkoqp6UXd/+K5+EQAAAAAArJ77Umj/XneffFcHVdWPZha9f5Dkmd192bT/DUkuTHJ4VR3Z3efMzeyT5LTMgvyTu3vrtP/3knwhyQlV9aHu/vzczEGZRfarkxzY3d+d9p+aZEuS06rqI0vnmhyZWWTfnOQ53X3bNHNmks8mOauqLuzu7+/g7wYAAAAAgJ10n3h0zA46PMnDkpyzFNmTZIraJ01vX70wc3SS3ZKcMR/Gp3j+1untqxZmlt6/ZSmyTzNbk7xrOt/LFmaWrnvSUmSfZr6Q2Z33D5vWDwAAAADA3eS+FNp3q6p/W1UnVtVvVdWztvG89WdP248v89lFSW5JclBV7bbCmY8tHLNTM1W1e5KDputfvAPXAQAAAABgDd2XHh3z40nev7Dv2qp6WXd/Zm7fY6ftVYsn6O47quraJD+dZN8kf7uCmW9W1c1JHlVVD+zuW6pqjySPTHJTd39zmbV+edruP7dvvyS7JLmmu+9Y4cw2VdWWbXz0uJXMAwAAAAAwc1+5o/29SZ6TWWzfI8kTkvy/kuyT5GNV9W/mjt1z2t6wjXMt7X/wTszsubBdi2s8eBufAwAAAACwBu4Td7R395sXdv1NkldV1U2ZfSHpyUl+9e5e13rq7gOW2z/d6b7pbl4OAAAAAMCGdV+5o31bzpy2B8/tW7z7fNHS/u/txMwNC9u1uMb3tvE5AAAAAABr4L4e2r81bfeY23fltP2hZ51X1a5JHp3kjiTXrHDm4dP5v97dtyRJd9+c5BtJHjR9vugx03b+me9XJ/lBkn2ndaxkBgAAAACANXZfD+1Pnbbz0fzCafv8ZY4/OMkDk2zu7ttXOPOChWN2aqa7b0uyebr+M3bgOgAAAAAArKF7fWivqsdX1R7L7N8nyRnT2w/MfXRekm8nObKqnjx3/O5JTpnevnvhdO9NcnuSY6bzLs08JMmJ09szF2aW3r9+Om5+Xb85ne+9CzNL1z1lWs/SzIFJjsjsDv0PLf6sAAAAAACsnfvCl6EekeSEqrooyVeTfD/Jfkl+McnuST6a5LSlg7v7xqp6RWbB/dNVdU6S65P8cpLHTvvPnb9Ad19bVa9L8o4kl1XVuUn+McnhSR6V5A+7+/MLM5ur6vQkxye5vKrOS3L/ab17JTm2u7cu/CznJHnhdN4vVdUFSR46zeyS5BXdfePO/qIAAAAAANhx94XQ/qnMAvnPJnl6Zs9L/16SzyZ5f5L3d3fPD3T3+VV1SJLXJ3lRZkH+K5lF8XcsHj/NvLOqtiZ5bZKXZvZ/C1yR5KTu/pPlFtbdJ1TVX2d2B/srk9yZ5ItJTu3ujyxzfFfVSzJ7hMzRSY5NcluSi5Kc0t2bV/5rAQAAAABgNdzrQ3t3fybJZ3Zi7nNJfmEHZy5IcsEOzpyd5OwdOP6OJG+bXgAAAAAArLN7/TPaAQAAAABgLQntAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAUI7AAAAAAAMENoBAAAAAGCA0A4AAAAAAAOEdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAAAAAABgjtAAAAAAAwQGgHAAAAAIABQjsAAAAAAAwQ2gEAAAAAYIDQDgAAAAAAA3Zd7wUAAAAAwD3FB/+//4/1XgI74H958aXrvQRI4o52AAAAAAAY4o52AAAA4D7n5JNPXu8lsAP89wLu6dzRDgAAAAAAA4R2AAAAAAAYILQDAAAAAMAAoR0AAAAAAAYI7QAAAAAAMEBoBwAAAACAAbuu9wIAAAAAAO7p/s15n1jvJbAD/urwn79br+eOdgAAAAAAGCC0AwAAAADAAKEdAAAAAAAGCO0AAAAAADBAaAcAAAAAgAFCOwAAAAAADBDaAQAAAABggNAOAAAAAAADhHYAAAAAABggtAMAAAAAwAChHQAA+L/bu/Mwyary8OPfl31RVkGQiYyCwKBRIsqqOEBEURSNmpiIMv7EaFwQDQpqBNwBQ0DEHXEUNAgqWyLixohsog4aVxZhVAQEZAbZZnDg/P54TzGXS1Vv1d3VPf39PE891X3udurWW3d577nnSpIkSeqDiXZJkiRJkiRJkvpgol2SJEmSJEmSpD6YaJckSZIkSZIkqQ8m2iVJkiRJkiRJ6oOJdkmSJEmSJEmS+mCiXZIkSZIkSZKkPpholyRJkiRJkiSpDybaJUmSJEmSJEnqg4l2SZIkSZIkSZL6YKJdkiRJkiRJkqQ+mGiXJEmSJEmSJKkPJtolSZIkSZIkSeqDiXZJkiRJkiRJkvpgol2SJEmSJEmSpD6YaJckSZIkSZIkqQ8m2iVJkiRJkiRJ6oOJdkmSJEmSJEmS+mCiXZIkSZIkSZKkPpholyRJkiRJkiSpDybaJUmSJEmSJEnqg4l2SZIkSZIkSZL6YKJdkiRJkiRJkqQ+mGiXJEmSJEmSJKkPJtolSZIkSZIkSeqDiXZJkiRJkiRJkvpgol2SJEmSJEmSpD6YaJckSZIkSZIkqQ8m2iVJkiRJkiRJ6oOJdkmSJEmSJEmS+mCiXZIkSZIkSZKkPpholyRJkiRJkiSpDybaJUmSJEmSJEnqg4l2SZIkSZIkSZL6YKJdkiRJkiRJkqQ+mGiXJEmSJEmSJKkPJtolSZIkSZIkSeqDifZpKCJmRcQpEXFjRCyLiEURcUJEbDjoukmSJEmSJEnSTLPaoCug0YmIrYBLgU2Bc4DfADsBbwGeGxG7l1L+PMAqSpIkSZIkSdKMYov26ecTZJL94FLKi0oph5dS9gKOB7YFPjjQ2kmSJEmSJEnSDGOifRqprdn3ARYBH28NPhK4G3hlRKw7yVWTJEmSJEmSpBnLRPv0smd9/1Yp5YHmgFLKncAlwDrALpNdMUmSJEmSJEmaqaKUMug6aIQi4iPAocChpZTjugw/CXgj8IZSyieHmddPegx6ytprr73qnDlzRlW3X//RbuGnkzlbbDxpy7rvpl9N2rLUvzU2337SlnXVLVdN2rLUv2033XbSlnXX1VdP2rI0Ph6xzTaTtqxbb7hj0pal/m0ya/1JW9bNi66btGWpf5vNfvykLWvpzXdO2rLUv7U2e+SkLeumm26atGWpf5tvvvmkLWvx4t9M2rLUvw033G7SlvXrxX+ZtGWpf3M2XG9M0y1cuPDLpZRXjHY6H4Y6vXTOVHqdYXbKN+hjGfffe++9dyxcuHBRH/NYWXS21CvdHnbhn3436CpMdyttbHDTwkHXYLpbaWNj4Q3GRp9W2tgAYKHx0aeVNj7+cMugazDtrbSxcePtbjf6tNLGBjcOugLT3kobG14Y6dtKGxvXX+8+pU8rbWwsvH5yl2eifYYqpew46DpMdZ1W/64rtRkb6sXYUC/GhoZifKgXY0O9GBvqxdhQL8aGejE2xo99tE8vnRbrve7B7ZQvmfiqSJIkSZIkSZLARPt00+nQuFdHqE+o73ZuK0mSJEmSJEmTxET79HJhfd8nIh7y3UXEI4HdgXuAyye7YpIkSZIkSZI0U5lon0ZKKb8FvgXMBt7YGvxeYF3g1FLK3ZNcNUmSJEmSJEmasXwY6vTzBuBS4MSI2Bv4NbAzsCfZZcy7B1g3SZIkSZIkSZpxopQy6DpolCLib4D3Ac8FNgZuAs4C3ltKWTzIukmSJEmSJEnSTGOiXZIkSZIkSZKkPthHuyRJkiRJkiRJfTDRLkmSJEmSJElSH0y0S5IkSZIkSZLUBxPtkiRJkiRJkiT1wUS7JEmSJEmSJEl9MNGugYqI2RFRImL+oOsy1UTEUXXdzB10XVZm03E91/ouGHQ9VnbGhnqZbrHhvnZyRcT8ur5nD7ouI2F8TI7ptt0A9ymTxdjQdIyB0XJfM/6m2/HGWETEvPoZ5w26LqMVEQsiogy6HjPNVFjvJtqlMYiIl0bExyLiBxHxl7rxP20E0+0WEd+IiNsj4t6I+L+IOCQiVp2Mek9HY1nXg1rPEbF7RBwbET+KiFsjYllEXB8RJ0fE1kNMt3ZEvDciroqIpRFxS0ScERFzJrK+01lEbBwRB0XEWRFxbf2e74iIiyPiNRHRdf82wNjYIyJOjYhfRMSf6/d8fUScGxF7DzGdsTEGEXFMRHw3Iv5Qv+fbI+LKiDgyIjbuMc2U2D5HxJo1TkpE3DDEeBtFxAkRsahua26MiFMiYtZk1ndlEBEH1PVdIuKgHuPsVw/c74iIuyLihxFx4CTU7ahG3bq9nttjOuNjlOq66rWeb+4xzUC3GxGxfkS8ry73rsjjpF9ExKcjYvUu47tP6UNE7F2PO25u/K4uiIjndRl3YLEREZtGHo/+IiLurMcdP4mIt0fEI3tMY2wMIyKeHxHfiogb6nd6XUScGRG79hh/UMecsyLi3bVu10bEA3U71vNcpE436hiYSfuamKRz0hjc8cac+v2fExG/b+z/VhtimgVD7DdLRKzVY7rta2zdUmPtqrrstSfuE2q8jfY3ESsutvV6nT7EtAdGxBX1N3FHjb39JuaT9a/nj0bSkP4DeApwF3ADsN1wE0TE/sDXgKXAV4DbgRcAxwO7Ay+bqMpOc6Na1wNez18DNgEuBb4ELAd2BV4DvDwinl1KuaxV3zWBb9e6/Rj4KPA3tZ7Pj4i9Sik/nMA6T1cvAz4J3ARcCPweeDTwD8DJwL4R8bJSyoNXswccG3vV1w+B7wF3A48FXgi8ICI+UEp5T3MCY6MvbwUWkuvvFmBdYBfgKOBfI2KXUsofOiNPse3zh4Athxoh8mLBpcA2ZDydTm4bX03Gxq6llOsmuqIrg4j4G+Akch/ziB7jvAn4GPBn4DTgPuClwPyI+NtSyqGTUNUvAIu6lF/bLjA++nIHcEKX8rvaBYPebkTEdsC3gC2A7wDnA6sDs4F/BP4d+GtjfPcpfYiIY4G3k8ei5wK3kcd8OwJzgW80xh1YbES2bP0hsCmwgIyLtYB9gGOBA+o+8N7GNMbGMCLiGOAd5H7gbPL73xrYH3hJRLyqlHJaY/xBbh+eBnwAKMD15HZtg6EmGEsMzMB9zYSfkw74eOM5wBHA/cA1td5dE+VdvLdH+fJ2QUTsTMbL6sBXgT+Q50hHAHtHxN6llGWjq7oGZNQ5sepn5Ha07RfdRo6I/ySPaW4APgusAbwcOC8i3lxKOWl01Z4EpRRfvgb2Ik8GCjB/0HUZZb33BJ4ABHlwXYDThhh/PTLZswx4WqN8LfIApQAvb01zVC2fO+jPO13W9aDXM3AY8Jgu5e+qy/h5l2HvrMPOBFZplO9fy3/ZLK/DCrBg0N/NgONiL/Jgtb1uNiOT7gV4yRSKjbV6lG8B/Ik8qN3c2Bi3+Oi1vj9Y19EnpkpstOY7F3gAeH2d/w09xvt0HX5cq/zgWv7NVvlspuG+dhLiJMgE5W+Bj9R1dFCXdbeUPOmd3SjfkExyF2DX1jTza/nscajjqGPN+Bjzul4ELBrhuIPep6wDXA0sBnbpMnw1IFpl7lPGvr5f2/mNAGt0Gb76FIqNj9d5HdkqXxX4bh32KmNjVOt0M/I47WZg09awPeu6uG4KxcAs4JnAevX/BXXeWw8xzVhiYEbta5j4c9LZDPZ4Y1tgZ2Dt+v+iOu/VhphmAVBGsYxVgV/V+b6wUb4KmXQvwOGtaebV8nmDjoExrNNRrZ/p9hrNb6KOP+ptALBbneZaYMPWvP5cfzOzp9p6t+sYTVkRsU5EvDMifhoRd9fbRC6LiH/uMm7U20kujeyyY2lklwEXRMQ/tcZ9ckT8d6y4xe3WiFgYedvbw26z7aaUcmEp5ZpSf8kj8FKy1cvppZQfN+azlLwSCPBvI5lRRDw2In4ZEfdFxCtHuPxpa5TreqDruZRyTCnlxi6DjgHuBZ4UjW4rIiLIhBrAO0opDzTmdQ7wA2B74FkjrPPbI28PvSQiNhrJNNNVKeV7pZTzmuuslt8MfKr+O7cxaNCxsbRH+R/JA+5VgMc3lmFs9KHX+gbOqO9PaJRNie1zRKxHnjB9t5TyqSHGewTwSvKuiKNag08Cfgc8JyIezzAiYpWI+Gi9XfPrMfNu2T2YvGj3anJ9dvP/gDWBk0opizqFpZTF5N0HsOK3OqSIeEpE/LHeXvvsMdd66GUYH5Nj0NuN15PbsXeWUi5vDyylLG8eN7lPGbva0veD5EX8fy2l3Ncep5Ty18a/g46Nzm/73FYd7wf+t/67SWMZxsbwtiSP035YSrmlOaCUciFwJ411yuCPOW8opfyglPKXES5j1DEwE/c1k3BOOtDjjVLKVaWUH5bG3S4T4FnAHOCiUsqD26gac++o/76+xuSQImLDiLiobm/eOTHVfdgy50XE1yK7jbq3rt9LIuKAYaZbMyI+ENlt6LKI+G1kd5Zr9Bh/74j4ZmSXQ8si4uqIODoi1m+N95u6bXhUj/kcVn9Xb2qVz4qIk+rnWBbZvdi5EfH00ayPMeTExqIT8x+sv4XOsheRF5bXJI/jhxURe0V2O3NjROww3hVtMtGuKSkiNgAuJncq9wOnkLdNbwJ8OSI+0Jrkg2SCYjMykfJfZCu1LWjclhURTyZvp9wfuLyOdwZwK/AG8oc6Efaq79/sMuwi4B5gt3ow31NEPAW4jLyV73mllFPHtZbT31Rdz4UVt87d3yjfiuw+5OpSyvVdpju/vu/VZdiD6gHqieQtwWcBe5dSbu+vytNa54S3ebvilIyNiNiUbD2yDLiqMcjYmBgvqO//1yibKrFxItlq6TXDjLcLsDZwSSnlzuaAeqJyQf13z2HquxbZcu1g8kD1pRN8cjWlRPY5ezTw0VLKRUOMOlR8jOh3WJe3NxlPAexRSvn2KKoL8IyIOLSeNP1Tr5MqjI9+rRnZZ/+7IuItEbFndO9Ld9DbjX8hjy1Oj+zz9N8iG6e8Iro/h8J9ytg9mzz/+DrwQGQ/3YfV+OjWN/egY+OX9f35rfmtAuxL3jX1vcYgY2N415BdeOzU3vZGxB7AI8nzzo5Bx8BojSUG3NcMbSwxMFWON0atHpccHhFvi4h9h4jtnp+xZDdDV5MXtoa8QBMRjwUuIePwVaWUD4+99qPySbJ+F5HdzJ1e/z81It4/xHRnkBdSziMvRBXyAtXX2hcVIuJ1rOjG6Wyyq6HbyTvnL615so4vkF3wPKwhanUgue36cmP+TwV+Sua+riK7KjoP2AO4OLo8c2QCPCYiXlePtV5Xc3S9jNfv4hV1/BvJO0N+OpoKj5Z9tGuqOgH4O+CwUsqxncK6Mz4beFdEfLXxA3kd8EfgSaWUe5ozah0QHUjesvWieoW+Od6G5E5vImxb369uDyilLI+I64EnkjuVX3ebQUT8PdnP293AM0spP5uguk5nU3U9v4w8CL+8lLJkJPWtrqnv2/Sacf1NfInsm/wk4C3tVt4zSeQDe15V/23ukKdEbETE04D9yP3vLDLxuz7w5lLKbSOpb2VsjEBEHEr2u70+2WfpM8gk+9GN0QYeGxHxYnL/dFAp5ffDjD4esbER2dpxN/IW3WNGUd1pr24nTiVbqL5rmNGHio+bIuJuYFZErNM+/mgs7wCywcC1wL6llN+NodrtE7hlEfER4IhWSyLjoz+bkbHRdH1EvLqU8v1G2cC2G5F3Xz6FbCTyWrJRSvOc7u6IOLiUcspI6lu5T+mt08JvKXAl8KTmwIi4iEwO3lqLBr1POZY8znh/ROxJPq9kDbKP9s3I/cyVjfGNjWGUUm6PiMPIBlq/ioizyS4LtiKftfNt8ly0Y9AxMFpjiQH3NUMbSwxMleONsWg/xPKWiHhjKeWrrfKRxM029fXbbiPUC1Dnk89fel4p5TvdxpsgTyqlPKReka3SzwcOj4hPlbxbuW0O8MROi+yIeDf5jLH9gAOoxx0RsSXZ8OYuYKdSym8ay/kEeRfEscC/1uJTyecxHEgmzJv1enpd7tc7Fz3r8e8Z5LnRns3jmoh4DPAj4HMRMbtMbD/5z66vZn0XAAc2z4MiYl2y4exdpZSbusxn2O1Mnc9hwIfJizP7T8ZFYFu0a8qpLXEOAH7cTLLDg7dbHUZeof2X1qR/5aGthTvT3NYuI7vxaI+3eAIPCju3+dzRY3infINuA+tO8xvkxYRdTLL3NOXWc0Q8jtzxLQfe1hrcb303IlvQvJi8KPXmlenEZoyOJk+Cv1FKuaBRPlVi42nAkcC7yYOi1YFXl1I+2RrP2Bgfh5Lr+xAyyf5NYJ9GQgQGHBsR8WjgM8D5pZTPjWCSfuu7JXmguRPwypXsxHakjiAv5s8bQcu5ka7v9bsNjIjDgS+Sd9PtPoaT3p+RraAeT7Ye3JJMri4hbz//4Bjru0GP+s7k+Pg8sDeZjFwX+FuyD+LZwPn15L5jkNuNjcjE+sbkieP7yRavjwIOIlvKnRwRzRZe7lPGbtP6/nZy3T6TbDzxZPJhtHuQLXI7BrpPKdm1yS5kq/K9yP3gwWSC6wwe2vJ6POo7I2KjlHICeSFhNXIbfDjZkOYPZH/DzS5lpsox50iNpb7ua4Y2ket0oo43xuIcstHQLPIYZTtyv7QB8JWIeG5r/H7j5tlkV0aFbK0/mUl22kn2WnYfeTfGauQxRDfvLw/t9mQp+VwEyGO8jgPIC6MnNZPs1bvJbqpe2bljoJRyA/nsjR0j4omt8Q+s719olD2fvED4sVbjAUp2f3sseQzU63P06x7ymGVH8i7eDcnuhC4ku3z9bk2ud/QbL6tExElkfuAs4NmTdaeVLdo1FT2dfFBGiYijugzv9KM+p1H2JeDNZCuDM4DvA5eVUto/yq8AbwHOjoivkgeGl3TbaE4hbyG7urmEfGjI4mHG19gMu54jYh55wt20oJSyoNsMI7sFOZ+85fiNpZTLxrG+j651fTxwQCnly8OMv9KLiIPJJ5L/huw3cryMW2yU7Hv7U7X11+PIfue+GBG7l1JG1O/iCBgbVSllM3gwmb0beaB1ZUTsV0pZOA6LGI/Y+Cx5PHbQONRnONuSt6GvS7Z0+u4kLHNKiYidyVbsx43zNrmb44EXka0SDyhdnh0QEYfw8BOEszt37JVSzmoN+z2ZRF1IdoF3aET8V49GBaM1o+OjlPLeVtEvyL5i7yL3LUeRycR+9bvd6DSUWhX4dCnlfY1xPhcR65At4g7joV2EjNVM36d01vdy8vtaVP//eb0b6SrgWRGx6zhsU/rep0TEbLKV8NrA8+q81qnzPQ7Yv9a1WxchozVjYiMi3kHePXIi2Wr/ZlYkFb8UETuUUt4xxCxGatzPR6agGb2vmSB9H2+MRSnl+FbRVWTPAzeSDc0+TPcuP8bipeSdOdeQcTPcHaDjLrLLmsPIRPRjye1s0xY9Jv1+l7KLyUaif9coe2p9f9i+u5SyOCKuJC/ubkc2xIDsPvnZZGL9HbWea5DdydxCXrTr6HR3tmWPPFvnGVZzWtONi3pB8ohW8UURsQ+5PnYmz4c+Ok6L/Br5u/gYcMhkXgQ20a6pqNO35NNZcbtmN49o/P1W4DryQQiH19fyiPgG8O+llGsBSilXRMQzySuCL6Um4yLiKuC9pZT/Hs8P0jDkFehG+ZIuw/YgW/B/1yT7sCZ6Pc+j+8OgFrQLapL9e+TB5FtKKZ8Y5/puRj7R/gZyxzSjRT7k5aPkk+y79Qs6ZWIDHmzJ8GvgLbVVwusi4juNWyyNjXFUSvkTcFZNTl5Ntvjp3P4/sNiIiFeRLYEOLN0fpNxNP/XdhmwN+1OyO4EZpd4y+0UyBt4zwsnuIFsKr092FdA2VGubPer7/3Q76a0OIVupNy0iv6OeSikLI+IKsg/PXcn+NZv1MD7Gz6fIRPsejbJB7lOasda+ENMpO5FsKdrhPmXsltT3KxtJdgBKKfdExAXkszV2IpOHgz7emE/ejfGUUkrneSR/AT5dL/KfQN7pNa8OMzaGERFzgWOAs0opzbtTF9aLLVcD/167jbiOwcfAaI2lvu5rhjbWdTrljjfG6GQy+b9DRDyyrOjHv5+42ZVscPlD8k6SSRX5YN8ryFbYPyDvaLqDTJbPJhPdvfqm/1O7oGQXQrex4q4pWPH5u3WT0izfoFF2FrmNPyAi3lnywdf7kb+xE0opzWeWdfJsL2Nojxhm+Liq6+JkMtG+BysS7f3EC3Vey4HzJvtOK7uO0VTU+UEdX0qJIV57diYopdxfSjmhlPIUsnXFS8iNzguBb0bjgRyllMtKKfuRG8ndydtXHk0+ZPXvJ+gzdR5y+LD+o+qJ/+PIjcB1XaZ9DblDOTIi3tdluFaY0PVcSpnbJQ6P6rKszcmD3e3Jluwnjra+Veeqcrd+7H5G7tC3IK8ED/nQmJVZbaHxMbL14Z6llJu7jDYlYqOHzoNc5o6kvpWxMQb1NtpfAU+MFc/vGGRsdFqufCEiSvNVy7dolG0wXH2roWLjPLI19w7k7ZndHpq4MnsEud7mAEtb6/vIOs5na9kJ9f+h4mNzskXeDaV7f6kvIvtJ/VxEvLZbhUops7vEx/wRfp5OF0jN22yNj/E3qvU80duNGmudJMOSLpN3knPNlnbuU8aus+6W9BjeXt8Di42IeCSZgL29kWRvurC+79goMzaGt199v7A9oP4eryDzKp2WqVP5mLObscSA+5qhjSUGpvLxxqjUZH8nuT5exyjvIu/WeTVwSuQDnifT28hE9Wvqb/DgUsp76m/vgqEn5dHtghoHjyKT5B2dPNhmPeazeWs8SnaBeEYd1un3vFu3Mc3p9h8mz9a+w28yPOxYq5RyN9mF1iPqb6BtqHiBfBjzYuDcmJyHvD7IRLumoiuAB8g+EEetlHJLKeXrpZR/JFsUb0XrwUV1vGWllEtLKUeQfRdC3qo3ETq3/7T7KYO80rYOcGnp/tCJJeRG8wfAeyLi2C7jKA18PUfELPL2sO2A1/doyd7xW7IbgG0i+3Jv27e+d731u5RyGvBy4DHkCc6QDwJZGUU+3OR4sjXGnuWhfWQ2DTw2htC5zbDZ4sDYmDiPqe+dZ3oMMjYuAz7X4wXZl2Hn/87yLyefM7J7Tao8qJ507FP/fVhCAKCU8mHyLrC/I1vVP+zgfyW2jN7r+8o6zsX1/04XEEPFx5C/QzIZugd5YvnpiHhjP5VvinwgZudCTfNE3fgYf7vU9+Z6HvQ+pdMv7cOObxtlza5B3KeM3XfJ/oC375HYaa/vQcbGGvV9vdp1QNsm9f2+RpmxMbxOg61Negxvr9dBbx9Gaywx4L5maGOJgSl5vDEWEbEt2ajxTqDZtV3Pz1gv1G0D/I7uF6GWkT0SnEne1XFaTVZPlq3r+9e6DOt2h8lww59BdgF3ZaOs8/fc9si1wc0O5IO52w9Rnl/fD4yITch4+b/y8G6BLq/vY8qzTbBux1rQx++iXnB+FplsPysiXtRnHUeulOLL18Be5G02hXyITLP8i7X8PcCqXabbCnhc/XtN8oEf7XFWJzdWBZhTy3YD1u4y7qF1vGPG8Bnm1mlPG2Kc9cirdMuApzXK1wIurdO/vDXNUbV8bv1/HfLEqgAfHfR3N6B4GXJdD3o9k7fjXUcm8OaNcJp31mWdCazSKN+/lv+yWV6HFbIvxs7/LyR3ujeRTzQf+Hc1SfHwnroufgxsNMy4g46NnXqUb0Xecl3IB7QYG/3HxTbA+l3KVyEfHFnIZ3NMidgY4nMUsuVSt2GfrsOPa5UfXMu/2SqfTWtfSz4f4AHymQaPGfT3NuhX4zs9qFX+uPob+jMwu1G+Idl6rAC7tqaZX8tn1/83IVt+FrI7u5HW6ZHAtl3K1yAfvFXIk632dsD4GP33PwdYt0v5bLI/2AK8q1E+6H3KjuSxxjXAJq3ld+Z5RGsa9yljj49z6rp4a6t8n/o7WUzd70yB2PhVne79rfK1yKRnAY41Nka1Tv+xft6bgS1aw/atMXAvsPFUiIEu9V9Q57P1EOOMJQZm7L6GiTknHdjxRo/PsKjOZ7Uewx9Hl/OvWofOZ/xMa9iqrNhGvbBRvkqNvQIc3ppmXi2f15jHqbXsa8Dqk/Sdf6ou8wWt8ueQjaUKcFRrWOe3dzWwYSsOLqvDXtX6jdxHXnDbujWvj9XxP9ujfleTDXQ658dv7TLO6jWW7gGe12M+uwLrjPdvoo7zVFrbkVq+d439AuzWGrZbLb+2tQ5n19/K0ubvpbneG/9vTV5M/CvwT5MSL5OxEF++er3onWhfr7HxuRo4hXyYxhfIFu8P7pzIPqoKebJxOtmHXqev5gKc05jv2eTtOf9LPsjmaPL2teXA7cBWI6z3i8gd23zyAR+FbA3QKfvPHtMsB+4i+y07ljyo6BzURGv8o2gccNWytWrdC3lwEyOp73R+jXZdD3I9k62ZCpn4ParHa3ZrmjXJBx4V4Ec1Jr9M7gjuBnbuspyHnNzUsueQO81byX45B/7dTXBcHFjXw3KyRXu3dT1vCsXGEnJb9hXgP8k+Us+r33MBTuwyjbExttg4hDzh/TbwGXLfcQq53ShkEmD7qRIbQ3yOQu9E+8Zkq6VCtrb8MLl/K2Q/kFu1xp9N933tPDJZdy3w2EF/dwOOm853elCXYW+uw24jE9zHk63HCt339/NpnPjWso3q77gA7x5hnWaTCYgryOOfo8nW9tfV+dwK7GB8jNv3f2f9DX+CPJb8at2WlFq+RmuagW43yAeKdb7Tz5In4VfXskuAtVrju08Ze3zMIk/SC5n8/EiNj+V1/b1kqsQG8Pdkcq+QrRf/C/gkK5Jm11ATwsbGiL//VchjikKeR36hbiPOJbfRhXwe05SIgTrt/MbrZlYkJTtlzxiHGJhR+xom+Jy0TjOQ44063aNacXNXnccXGmXbtb63peQ28TONmFnSiKMNuixn5xpT99Xxj27U92JgzS7xUWic15G/yc/W8vPa00zQ9/9kctu6FDitfp/fILcBpzN0ov0c4Eby+SnHseLCyf+04wB4Ayu2NSeTv6vOhYtf06NxGfAfdZz7yN/tpkN8jptYcazwcXKfdjorzpU2m6DfxAKyK5gzydg+nhV3jRXgP3os57g6/A91mo+Tv5ECvKnL+AtoJNpr2Za1bstpXNyYsHiZ6AX48jXUix473DpsDeBNdcNyR92w/b7+GA9hRauB1cknLJ9fhy8lD/AuJ6+Sr9GY5z7A58kk/B3kRv6qutHbchT1PqqxQej2WtRjut3JDfJi8uTt5+Rtc91a7XeWMbfLevl6Z73R5argyvQay7oe1Hoepp6l23LqdOsA7yNPfpbV+D2TVjKwtZwFXcrnkomC24GnD/q7G3Bc9FpHg4qNg8mDqd+RJ6Gd7dmZwHOGmM7YGH1sPIm8kPpT8iBsObm9/1H9TnsdoE6p7TNDJNrr8I3Ii8q/Iw+qbyIvKMzqMu7sTp26DPtn8oB8EfD4QX9/A4ybznf6sER7Hf4CsluwO8ljhx+RD7HtNu58Wie+tXw9ViQy3j+COq1HHp9cTiZK7iNPfH9Gnph2PYkyPsb0/T8L+G8yCbKkfuZbyeTaq+iR4Br0dgP4B+Ai8oR8Kdnq9N30SDrgPqWfGNmEvJjR+U3dRj4PqtcdawOLDTKRcip5nHFfXf4vgQ/RJfFlbIxona5Onn9eXn9vy4FbyGO7faZgDAx3jDyv3xio08yYfQ0TfE7amGbSjzda38+IzmPJhy7Pr5/pz/X7u53s9ujNtC5Ot5a1fY2t22qsXQ28l+49D8zrFrPkQ4M7rbwv6DbtBMTAbmQ3JYvr93MxmWyey9CJ9jWBD5CN8paRDSaOpPe+eh/yYauL6/jXkon9DYao22PJC1aFfPjnUJ9jU/I48hfkeeld5O/+q8AB9LiLod/fBPkMiv8hf+d3seJ8+CvAM4dZ1rz6W7i7rvvvA/v1GHcBrUR7Ld+CPM67H3jtRMZK1AVKkiRJkiRJkqQx8GGokiRJkiRJkiT1wUS7JEmSJEmSJEl9MNEuSZIkSZIkSVIfTLRLkiRJkiRJktQHE+2SJEmSJEmSJPXBRLskSZIkSZIkSX0w0S5JkiRJkiRJUh9MtEuSJEmSJEmS1AcT7ZIkSZIkSZIk9cFEuyRJkiRJkiRJfTDRLkmSJEmSJElSH0y0S5IkSQIgImZHRImI+YOuiyRJkjSdmGiXJEmSNKNFxLx6gWHeoOsiSZKk6Wm1QVdAkiRJ0pTxR2AOcMegKyJJkiRNJybaJUmSJAFQSvkr8JtB10OSJEmabuw6RpIkSRLQvY/2iJhfy2ZHxOsi4ucRsTQi/hQRn4mI9XvMa1ZEnBgR10TEvRFxe0RcERHv6aN+G0XEByPiFxFxT0TcERE/i4ijI2Ldxng7RsRH67Dba32viYjjImLD1jwXAJ+v/36+ftbOa/ZY6ypJkqSZxRbtkiRJkkbiWOA5wHnAt4A9gdcCWwN7NUeMiKcBFwAbARcBXwfWAbYHjgLeP9qFR8TjgAuBLYGfAJ8kGw5tA7wV+BRwdx39tcCLge8D36nj7Qi8Ddg3InYupdxZx50PLAH2B84BftpY7JLR1lOSJEkzk4l2SZIkSSOxC/C3pZTfA0TEasD3gD0jYqdSyhW1fA3gTDLJ/opSypebM4mIWWNc/pfIJPu7Sikfbs3zUcBdjaIPA28spdzfGu81wMnAG4BjAEop8yMCMtF+dill/hjrJ0mSpBnMrmMkSZIkjcT7Okl2gFLKclZ0ubJTY7wXALOBc9tJ9jrdDaNdcETsCOxKtjY/pss8byulLG38/7t2kr06BfgL2TJfkiRJGjcm2iVJkiSNxI+7lP2hvjf7Pd+lvp8/jsvuzPOCUsoDw40cEatHxJsi4uLaR/v9EVGAB4D1gC3GsW6SJEmSXcdIkiRJGpElXcqW1/dVG2Ub1Pc/juOyRzvPr5B9tF9H9rt+M7CsDjsEWHMc6yZJkiSZaJckSZI0rpbU9/FsNT7iedYHsb6YfAjqvrWLm86wVYB3jGO9JEmSJMCuYyRJkiSNr8vr+74TMM/n1GT5ULau7+c2k+zVTsDaXabp9Oe+apdhkiRJ0rBMtEuSJEkaT+cBi4AXRsQ/twdGxKzRzrCU8hPgUmAH4LAu89w4Itaq/y6q73Nb42wKfLzHIv5c3x872rpJkiRJYNcxkiRJksZRKeW+iHgZ8C3gyxHxOrJF+lrAHGBvxnYecgCwAPhQRLyk/h3AE4B9gO3IJPuPgEuAf4iIS4GLgUeTLeyvAm7sMu/LgHuAQyJiY7JPd4CPlVLuGENdJUmSNMOYaJckSZI0rkopP46IHYDDyQT3bsCdwLXAEWOc5/UR8VSyj/UXAW8ClpLJ9eOAW+p490fEC4EPAM8DDiYfonpyLftVl3kvrsn7I4F5wLp10GmAiXZJkiQNK0opg66DJEmSJEmSJEnTln20S5IkSZIkSZLUBxPtkiRJkiRJkiT1wT7aJUmSJE26iDgE2GAEoy4opSyY0MpIkiRJfTLRLkmSJGkQDgG2HOG4CyauGpIkSVL/fBiqJEmSJEmSJEl9sI92SZIkSZIkSZL6YKJdkiRJkiRJkqQ+mGiXJEmSJEmSJKkPJtolSZIkSZIkSeqDiXZJkiRJkiRJkvpgol2SJEmSJEmSpD6YaJckSZIkSZIkqQ8m2iVJkiRJkiRJ6oOJdkmSJEmSJEmS+mCiXZIkSZIkSZKkPpholyRJkiRJkiSpDybaJUmSJEmSJEnqg4l2SZIkSZIkSZL68P8BQPvJRLeh5RwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 756x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 495,
       "width": 749
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = sns.catplot(x=\"inc_cat\", kind=\"count\", data=df, height=7, aspect=1.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Draw a sample for conducting the DP analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 10000\n",
    "\n",
    "# draw sample \n",
    "sample = df.sample(n = sample_size, random_state=0)\n",
    "# save sample for DP analysis\n",
    "col_names = [col for col in sample.columns]\n",
    "sample.to_csv(samplefile, index=False, header=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define two functions to support the further analysis:\n",
    "1. `dp_histo` is used to generate data for the differentially private histograms for a given privacy parameter epsilon. It generates counts based on two alternative mechanisms to generate statistical noise: Geometric and Laplace.\n",
    "2. `perc_error` is used to calculate the percentage deviation between the true and the differentially private distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "ename": "OpenDPException",
     "evalue": "MetricMismatch(\"Intermediate metrics don't match. See https://github.com/opendp/opendp/discussions/297\n    output_metric: L1Distance()\n    input_metric:  SymmetricDistance()\n\")\n\t   0: backtrace::backtrace::trace\n\t   1: backtrace::capture::Backtrace::new_unresolved\n\t   2: _opendp_comb__make_chain_tt\n\t   3: _ffi_call_unix64\n\t   4: _ffi_call_int\n\t   5: __ctypes_callproc\n\t   6: _PyCFuncPtr_call\n\t   7: __PyObject_MakeTpCall\n\t   8: _call_function\n\t   9: __PyEval_EvalFrameDefault\n\t  10: __PyFunction_Vectorcall\n\t  11: _call_function\n\t  12: __PyEval_EvalFrameDefault\n\t  13: __PyFunction_Vectorcall\n\t  14: _slot_nb_rshift\n\t  15: _binary_op1\n\t  16: __PyEval_EvalFrameDefault\n\t  17: __PyEval_EvalCodeWithName\n\t  18: _builtin_exec\n\t  19: _cfunction_vectorcall_FASTCALL\n\t  20: _call_function\n\t  21: __PyEval_EvalFrameDefault\n\t  22: _gen_send_ex\n\t  23: __PyEval_EvalFrameDefault\n\t  24: _gen_send_ex\n\t  25: __PyEval_EvalFrameDefault\n\t  26: _gen_send_ex\n\t  27: _method_vectorcall_O\n\t  28: _call_function\n\t  29: __PyEval_EvalFrameDefault\n\t  30: __PyFunction_Vectorcall\n\t  31: _call_function\n\t  32: __PyEval_EvalFrameDefault\n\t  33: __PyFunction_Vectorcall\n\t  34: _call_function\n\t  35: __PyEval_EvalFrameDefault\n\t  36: __PyEval_EvalCodeWithName\n\t  37: __PyFunction_Vectorcall\n\t  38: _method_vectorcall\n\t  39: _PyVectorcall_Call\n\t  40: __PyEval_EvalFrameDefault\n\t  41: __PyEval_EvalCodeWithName\n\t  42: __PyFunction_Vectorcall\n\t  43: _method_vectorcall\n\t  44: _call_function\n\t  45: __PyEval_EvalFrameDefault\n\t  46: _gen_send_ex\n\t  47: __PyEval_EvalFrameDefault\n\t  48: _gen_send_ex\n\t  49: __PyEval_EvalFrameDefault\n\t  50: _gen_send_ex\n\t  51: __PyEval_EvalFrameDefault\n\t  52: _gen_send_ex\n\t  53: __PyEval_EvalFrameDefault\n\t  54: _gen_send_ex\n\t  55: _task_step\n\t  56: _TaskWakeupMethWrapper_call\n\t  57: __PyObject_MakeTpCall\n\t  58: _context_run\n\t  59: _cfunction_vectorcall_FASTCALL_KEYWORDS\n\t  60: _PyVectorcall_Call\n\t  61: __PyEval_EvalFrameDefault\n\t  62: __PyFunction_Vectorcall\n\t  63: _call_function\n\t  64: __PyEval_EvalFrameDefault\n\t  65: __PyFunction_Vectorcall\n\t  66: _call_function\n\t  67: __PyEval_EvalFrameDefault\n\t  68: __PyFunction_Vectorcall\n\t  69: _call_function\n\t  70: __PyEval_EvalFrameDefault\n\t  71: __PyFunction_Vectorcall\n\t  72: _call_function\n\t  73: __PyEval_EvalFrameDefault\n\t  74: __PyFunction_Vectorcall\n\t  75: _call_function\n\t  76: __PyEval_EvalFrameDefault\n\t  77: __PyEval_EvalCodeWithName\n\t  78: __PyFunction_Vectorcall\n\t  79: _method_vectorcall\n\t  80: _call_function\n\t  81: __PyEval_EvalFrameDefault\n\t  82: __PyEval_EvalCodeWithName\n\t  83: _builtin_exec\n\t  84: _cfunction_vectorcall_FASTCALL\n\t  85: _call_function\n\t  86: __PyEval_EvalFrameDefault\n\t  87: __PyEval_EvalCodeWithName\n\t  88: __PyFunction_Vectorcall\n\t  89: _call_function\n\t  90: __PyEval_EvalFrameDefault\n\t  91: __PyEval_EvalCodeWithName\n\t  92: __PyFunction_Vectorcall\n\t  93: _PyVectorcall_Call\n\t  94: _pymain_run_module\n\t  95: _pymain_run_python\n\t  96: _Py_RunMain\n\t  97: _pymain_main\n\t  98: _main\n\t",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mOpenDPException\u001b[0m                           Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/nx/9gd39tcx2lq7c3vmpmfk8vjc0000gn/T/ipykernel_61437/2199303115.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     54\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     55\u001b[0m lap_histogram = (\n\u001b[0;32m---> 56\u001b[0;31m     \u001b[0mmake_split_dataframe\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseparator\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\",\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcol_names\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcol_names\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m>>\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     57\u001b[0m     \u001b[0mmake_select_column\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"inc_cat\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTOA\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m>>\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     58\u001b[0m     \u001b[0mmake_count_by_categories\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcategories\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mincome_categories\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m>>\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/envs/synth/lib/python3.8/site-packages/opendp/mod.py\u001b[0m in \u001b[0;36m__rshift__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    197\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTransformation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    198\u001b[0m             \u001b[0;32mfrom\u001b[0m \u001b[0mopendp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcomb\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmake_chain_tt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 199\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mmake_chain_tt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    200\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    201\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"rshift expected a measurement or transformation, got {other}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/envs/synth/lib/python3.8/site-packages/opendp/comb.py\u001b[0m in \u001b[0;36mmake_chain_tt\u001b[0;34m(transformation1, transformation0)\u001b[0m\n\u001b[1;32m     71\u001b[0m     \u001b[0mfunction\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrestype\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFfiResult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     72\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 73\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mc_to_py\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0munwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtransformation1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtransformation0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTransformation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     74\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     75\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/envs/synth/lib/python3.8/site-packages/opendp/_lib.py\u001b[0m in \u001b[0;36munwrap\u001b[0;34m(result, type_)\u001b[0m\n\u001b[1;32m    127\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mOpenDPException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Failed to free error.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    128\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 129\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mOpenDPException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvariant\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmessage\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbacktrace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mOpenDPException\u001b[0m: MetricMismatch(\"Intermediate metrics don't match. See https://github.com/opendp/opendp/discussions/297\n    output_metric: L1Distance()\n    input_metric:  SymmetricDistance()\n\")\n\t   0: backtrace::backtrace::trace\n\t   1: backtrace::capture::Backtrace::new_unresolved\n\t   2: _opendp_comb__make_chain_tt\n\t   3: _ffi_call_unix64\n\t   4: _ffi_call_int\n\t   5: __ctypes_callproc\n\t   6: _PyCFuncPtr_call\n\t   7: __PyObject_MakeTpCall\n\t   8: _call_function\n\t   9: __PyEval_EvalFrameDefault\n\t  10: __PyFunction_Vectorcall\n\t  11: _call_function\n\t  12: __PyEval_EvalFrameDefault\n\t  13: __PyFunction_Vectorcall\n\t  14: _slot_nb_rshift\n\t  15: _binary_op1\n\t  16: __PyEval_EvalFrameDefault\n\t  17: __PyEval_EvalCodeWithName\n\t  18: _builtin_exec\n\t  19: _cfunction_vectorcall_FASTCALL\n\t  20: _call_function\n\t  21: __PyEval_EvalFrameDefault\n\t  22: _gen_send_ex\n\t  23: __PyEval_EvalFrameDefault\n\t  24: _gen_send_ex\n\t  25: __PyEval_EvalFrameDefault\n\t  26: _gen_send_ex\n\t  27: _method_vectorcall_O\n\t  28: _call_function\n\t  29: __PyEval_EvalFrameDefault\n\t  30: __PyFunction_Vectorcall\n\t  31: _call_function\n\t  32: __PyEval_EvalFrameDefault\n\t  33: __PyFunction_Vectorcall\n\t  34: _call_function\n\t  35: __PyEval_EvalFrameDefault\n\t  36: __PyEval_EvalCodeWithName\n\t  37: __PyFunction_Vectorcall\n\t  38: _method_vectorcall\n\t  39: _PyVectorcall_Call\n\t  40: __PyEval_EvalFrameDefault\n\t  41: __PyEval_EvalCodeWithName\n\t  42: __PyFunction_Vectorcall\n\t  43: _method_vectorcall\n\t  44: _call_function\n\t  45: __PyEval_EvalFrameDefault\n\t  46: _gen_send_ex\n\t  47: __PyEval_EvalFrameDefault\n\t  48: _gen_send_ex\n\t  49: __PyEval_EvalFrameDefault\n\t  50: _gen_send_ex\n\t  51: __PyEval_EvalFrameDefault\n\t  52: _gen_send_ex\n\t  53: __PyEval_EvalFrameDefault\n\t  54: _gen_send_ex\n\t  55: _task_step\n\t  56: _TaskWakeupMethWrapper_call\n\t  57: __PyObject_MakeTpCall\n\t  58: _context_run\n\t  59: _cfunction_vectorcall_FASTCALL_KEYWORDS\n\t  60: _PyVectorcall_Call\n\t  61: __PyEval_EvalFrameDefault\n\t  62: __PyFunction_Vectorcall\n\t  63: _call_function\n\t  64: __PyEval_EvalFrameDefault\n\t  65: __PyFunction_Vectorcall\n\t  66: _call_function\n\t  67: __PyEval_EvalFrameDefault\n\t  68: __PyFunction_Vectorcall\n\t  69: _call_function\n\t  70: __PyEval_EvalFrameDefault\n\t  71: __PyFunction_Vectorcall\n\t  72: _call_function\n\t  73: __PyEval_EvalFrameDefault\n\t  74: __PyFunction_Vectorcall\n\t  75: _call_function\n\t  76: __PyEval_EvalFrameDefault\n\t  77: __PyEval_EvalCodeWithName\n\t  78: __PyFunction_Vectorcall\n\t  79: _method_vectorcall\n\t  80: _call_function\n\t  81: __PyEval_EvalFrameDefault\n\t  82: __PyEval_EvalCodeWithName\n\t  83: _builtin_exec\n\t  84: _cfunction_vectorcall_FASTCALL\n\t  85: _call_function\n\t  86: __PyEval_EvalFrameDefault\n\t  87: __PyEval_EvalCodeWithName\n\t  88: __PyFunction_Vectorcall\n\t  89: _call_function\n\t  90: __PyEval_EvalFrameDefault\n\t  91: __PyEval_EvalCodeWithName\n\t  92: __PyFunction_Vectorcall\n\t  93: _PyVectorcall_Call\n\t  94: _pymain_run_module\n\t  95: _pymain_run_python\n\t  96: _Py_RunMain\n\t  97: _pymain_main\n\t  98: _main\n\t"
     ]
    }
   ],
   "source": [
    "from opendp.meas import *\n",
    "from opendp.mod import enable_features, binary_search_chain, Measurement, Transformation\n",
    "from opendp.trans import *\n",
    "from opendp.typing import *\n",
    "\n",
    "# the greatest number of records that any one individual can influence in the dataset\n",
    "max_influence = 1\n",
    "# we can also reasonably intuit that age and income will be numeric,\n",
    "#     as well as bounds for them, without looking at the data\n",
    "age_bounds = (0, 100)\n",
    "income_bounds = (0, 150_000)\n",
    "\n",
    "budget = (0.1, 1e-8)\n",
    "\n",
    "from opendp.mod import enable_features\n",
    "enable_features('contrib', 'floating-point')\n",
    "\n",
    "from opendp.trans import make_split_dataframe, make_select_column\n",
    "\n",
    "histogram = (\n",
    "    # Convert data into a dataframe where columns are of type Vec<str>\n",
    "    make_split_dataframe(separator=\",\", col_names=col_names) >>\n",
    "    # Selects a column of df, Vec<str>\n",
    "    make_select_column(key=\"inc_cat\", TOA=str) >>\n",
    "    make_count_by_categories(categories=income_categories)\n",
    ")\n",
    "\n",
    "\n",
    "geom_histogram = binary_search_chain(\n",
    "    lambda s: histogram >> make_base_geometric(scale=s, D=VectorDomain[AllDomain[int]]),\n",
    "    d_in=max_influence, d_out=budget[0])\n",
    "\n",
    "#lap_histogram = binary_search_chain(\n",
    "#    lambda s: histogram >>\n",
    "#    make_cast_default(TIA=int, TOA=float) >>\n",
    "#    make_base_laplace(scale=1.0),\n",
    "#    d_in=max_influence, d_out=budget[0])\n",
    "\n",
    "#hl1 = make_split_dataframe(separator=\",\", col_names=col_names) >> make_select_column(key=\"inc_cat\", TOA=str)\n",
    "\n",
    "#hl2 = hl1 >> make_count_by_categories(categories=income_categories)\n",
    "\n",
    "#hl3 = hl2 >> make_cast_default(TIA=int, TOA=float) \n",
    "\n",
    "#hl4 = hl3 >> make_base_laplace(scale=1.0, D=VectorDomain[AllDomain[int]])\n",
    "\n",
    "\n",
    "geom_histogram = (\n",
    "    make_split_dataframe(separator=\",\", col_names=col_names) >>\n",
    "    make_select_column(key=\"inc_cat\", TOA=str) >>\n",
    "    make_count_by_categories(categories=income_categories) >>\n",
    "    make_base_geometric(scale=1.0, D=VectorDomain[AllDomain[int]])\n",
    ")\n",
    "\n",
    "lap_histogram = (\n",
    "    make_split_dataframe(separator=\",\", col_names=col_names) >>\n",
    "    make_select_column(key=\"inc_cat\", TOA=str) >>\n",
    "    make_count_by_categories(categories=income_categories) >>\n",
    "    make_cast(TIA=int, TOA=float) >>\n",
    "    make_base_laplace(scale=1.0, D=VectorDomain[AllDomain[float]])\n",
    ")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def perc_error(true_col, dp_col):\n",
    "    return sum(abs(true_col - dp_col)) / sum(true_col) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Income histogram actual counts:   [3410, 1849, 1316, 992, 656, 494, 518, 279, 257, 229, 0]\n",
      "Income histogram Geometric DP release:   [3420, 1859, 1316, 999, 658, 487, 518, 288, 262, 223, 19]\n"
     ]
    }
   ],
   "source": [
    "with open(samplefile) as input_data:\n",
    "    data = input_data.read()\n",
    "\n",
    "actual = histogram(data)\n",
    "geo_counts = geom_histogram(data) \n",
    "\n",
    "print(\"Income histogram actual counts:   \" + str(actual))\n",
    "print(\"Income histogram Geometric DP release:   \" + str(geo_counts))\n",
    "#print(\"Income histogram Laplace DP release:     \" + str(laplace_counts.astype(int)))\n",
    "\n",
    "# Clarify meaning of 11th value (bins = 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = pd.DataFrame(index=income_categories)\n",
    "private_df = pd.read_csv(samplefile)\n",
    "results['True'] = private_df['inc_cat'].value_counts()\n",
    "results['Geom-DP'] = geo_counts[:len(income_categories)]\n",
    "results['Lapl-DP'] = laplace_counts[:len(income_categories)].astype(int)\n",
    "\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_geom = perc_error(results['True'], results['Geom-DP'])\n",
    "error_lapl = perc_error(results['True'], results['Lapl-DP'])\n",
    "\n",
    "print ('Geometric DP mechanism with ε = ' + str(epsilon) + ' : Deviation compared to non-private sample  = ' + str(round(error_geom,1)) +'%')\n",
    "print ('Laplace DP mechanism with ε = ' + str(epsilon) + ' : Deviation compared to non-private sample = ' + str(round(error_lapl,1)) +'%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assessing the impact of privacy guarantee and amount of data on accuracy\n",
    "We will generate histograms for various sample sizes. Each histogram compares the data from the original distribution with differentially private counts at different privacy levels (controlled by the privacy parameter ε). Low values of ε are associated with higher privacy guarantees and therefore higher amounts of statistical noise. The following charts can be used to develop an intuition for the tradeoff between privacy and accuracy.\n",
    "Furthermore, we will see how the loss in accuracy can be compensated by increasing the sample size. Feel free to adjust `sample_sizes` and `epsilons` in the cell below to investigate different settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_sizes = [10000, 50000, 200000]\n",
    "epsilons = [0.5, 0.1, 0.05]\n",
    "\n",
    "for sample_size in sample_sizes:\n",
    "    # draw sample \n",
    "    sample = df.sample(n = sample_size, random_state=0)\n",
    "    # save sample for DP analysis\n",
    "    sample.to_csv(samplefile, index=False)\n",
    "\n",
    "    results = pd.DataFrame(index=income_categories)\n",
    "    private_df = pd.read_csv(samplefile)\n",
    "    results['True'] = private_df['inc_cat'].value_counts()\n",
    "\n",
    "    labels = ['True value from original sample']\n",
    "    \n",
    "    for epsilon in epsilons:\n",
    "        geo_counts, laplace_counts = dp_histo(epsilon)\n",
    "\n",
    "        results['Geom-DP-'+str(epsilon)] = geo_counts[:len(income_categories)]\n",
    "               \n",
    "        error = perc_error(results['True'], results['Geom-DP-'+str(epsilon)])\n",
    "            \n",
    "        labels.append('DP with ε = ' + str(epsilon) + ' : Error = ' + str(round(error,1)) + '%')\n",
    "\n",
    "    ax = results.plot.bar(rot=45, figsize=(13, 8), width=0.8, fontsize=12, colormap = 'Accent')\n",
    "    plt.title('Income Distribution (n = '+str(sample_size)+')', fontsize=16)\n",
    "    plt.legend(fontsize = 12, labels=labels)\n",
    "    plt.grid(axis='y', alpha=0.5)\n",
    "    plt.ylabel('Count', fontsize = 14)\n",
    "    plt.ylim(-50)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
